Set up libraries

library(org.Mm.eg.db)
library(AnnotationDbi)
library(data.table)
library(tidyverse)
library(rtracklayer)
library(RColorBrewer)
library(gplots)
library(clusterProfiler)
library(fgsea)
library(grid)

KEGG analysis

dir.create("PDF_figure/KEGG_and_GSEA_analysis_merged", showWarnings = FALSE)

# lead the miRNA list with peaks associated
mirs.peaks <- readRDS("Datafiles/miRNA-merged-peaks-list-09282019-withIDs.rds")
len <- sapply(mirs.peaks, function(x) length(x))
mirs.peaks <- mirs.peaks[order(-len)]
mirna <- names(mirs.peaks)

Run KEGG analysis with the whole universe as background.

KEGG.list <- list()
mirname <- c()
for (i in 1:20) {
  sig.genes <- mirs.peaks[[i]]$target_Entrez_ID
  kk <- enrichKEGG(gene = sig.genes, organism = 'mmu')
  enriched.number <- dim(kk)[1]
  print(paste("Number of KEGG pathway enriched for", mirna[i],  ":", enriched.number))
  if (enriched.number > 0) {
    KEGG.list <- c(KEGG.list, list(as.data.frame(kk)))
    d <- dotplot(kk, showCategory = 20) + labs(title = paste("Enriched KEGG pathways for", mirna[i], "targets"))
    print(d)
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(kk), paste("KEGG Analysis/", str_remove_all(mirna[i], "/"), "_KEGG.csv", sep = ""))
  }
}
## Reading KEGG annotation online:
## 
## Reading KEGG annotation online:
## [1] "Number of KEGG pathway enriched for let-7-5p/miR-98-5p : 1"

## [1] "Number of KEGG pathway enriched for miR-15-5p/16-5p/195-5p/322-5p/497-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-29-3p : 11"

## [1] "Number of KEGG pathway enriched for miR-200-3p/429-3p : 60"

## [1] "Number of KEGG pathway enriched for miR-17-5p/20-5p/93-5p/106-5p : 1"

## [1] "Number of KEGG pathway enriched for miR-24-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-194-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-27-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-103-3p/107-3p : 1"

## [1] "Number of KEGG pathway enriched for miR-196-5p : 1"

## [1] "Number of KEGG pathway enriched for miR-23-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-26-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-19-3p : 3"

## [1] "Number of KEGG pathway enriched for miR-25-3p/32-5p/92-3p/363-3p/367-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-30-5p/384-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-484 : 3"

## [1] "Number of KEGG pathway enriched for miR-141-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-130-3p/301-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-22-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-96-5p : 0"
pdf("PDF_figure/KEGG_and_GSEA_analysis_merged/KEGG.pdf",
    height = 6,
    width = 10)
for (i in 1:20) {
  sig.genes <- mirs.peaks[[i]]$target_Entrez_ID
  kk <- enrichKEGG(gene = sig.genes, organism = 'mmu')
  enriched.number <- dim(kk)[1]
  print(paste("Number of KEGG pathway enriched for", mirna[i],  ":", enriched.number))
  if (enriched.number > 0) {
    KEGG.list <- c(KEGG.list, list(as.data.frame(kk)))
    d <- dotplot(kk, showCategory = 20) + labs(title = paste("Enriched KEGG pathways for", mirna[i], "targets"))
    print(d)
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(kk), paste("KEGG Analysis/", str_remove_all(mirna[i], "/"), "_KEGG.csv", sep = ""))
  }
}
## [1] "Number of KEGG pathway enriched for let-7-5p/miR-98-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-15-5p/16-5p/195-5p/322-5p/497-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-29-3p : 11"
## [1] "Number of KEGG pathway enriched for miR-200-3p/429-3p : 60"
## [1] "Number of KEGG pathway enriched for miR-17-5p/20-5p/93-5p/106-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-24-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-194-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-27-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-103-3p/107-3p : 1"
## [1] "Number of KEGG pathway enriched for miR-196-5p : 1"
## [1] "Number of KEGG pathway enriched for miR-23-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-26-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-19-3p : 3"
## [1] "Number of KEGG pathway enriched for miR-25-3p/32-5p/92-3p/363-3p/367-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-30-5p/384-5p : 0"
## [1] "Number of KEGG pathway enriched for miR-484 : 3"
## [1] "Number of KEGG pathway enriched for miR-141-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-130-3p/301-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-22-3p : 0"
## [1] "Number of KEGG pathway enriched for miR-96-5p : 0"
dev.off()
## quartz_off_screen 
##                 2
names(KEGG.list) <- mirname
saveRDS(KEGG.list, "Datafiles/merged.peak.KEGG.result.rds")

GSEA analysis

I would like to do GSEA analysis on KEGG pathways for each miRNA's associated peaks.

Hallmark geneset

# load the proteomics datafile
protein_DGE <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/Proteomics data/colon tumor-enema model/crcMS_diff.csv")
## Warning: Missing column names filled in: 'X1' [1]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   `Protein Id` = col_character(),
##   Description = col_character(),
##   p_values = col_double(),
##   q_values = col_double(),
##   foldChange = col_double(),
##   LFC = col_double()
## )
protein_quant <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/Proteomics data/colon tumor-enema model/2017-03-19_Haigis_5v5_Pro.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Protein Id` = col_character(),
##   Description = col_character(),
##   control1 = col_double(),
##   control2 = col_double(),
##   control3 = col_double(),
##   control4 = col_double(),
##   control5 = col_double(),
##   KRAS1 = col_double(),
##   KRAS2 = col_double(),
##   KRAS3 = col_double(),
##   KRAS4 = col_double(),
##   KRAS5 = col_double()
## )
# annotate the proteomics dataset with Entrez ID
# Entrez IDs are annotated using `org.Mm.eg.db` package since this is the most updated
annotations_entrez <- AnnotationDbi::select(org.Mm.eg.db,
                                           keys = as.character(protein_DGE$`Protein Id`),
                                           columns = c("ENTREZID"),
                                           keytype = "UNIPROT")
## 'select()' returned 1:many mapping between keys and columns
# Determine the indices for the non-duplicated genes
non_duplicates_idx <- which(duplicated(annotations_entrez$ENTREZID) == FALSE)

# Return only the non-duplicated genes using indices
annotations_entrez <- annotations_entrez[non_duplicates_idx, ]

annotations_entrez <- annotations_entrez[!is.na(annotations_entrez$ENTREZID),]

# Check number of NAs returned
is.na(annotations_entrez$ENTREZID) %>%
  which() %>%
  length()
## [1] 0
# annotate the dataset with Entrez ID
protein_DGE$Entrez_ID <- NA
for (i in 1:dim(protein_DGE)[1]) {
  index <- grep(protein_DGE$`Protein Id`[i], annotations_entrez$UNIPROT)
  if (length(index)>0) {
    protein_DGE$Entrez_ID[i] <- paste(annotations_entrez$ENTREZID[index], collapse = " ")
  }
}

# calculate signal-to-noise ratio for GSEA later
s2n <- function(num_list, cond_1 = c(4:6), cond_2 = c(1:3)) {
  mean1 <- mean(as.numeric(num_list[cond_1]))
  if (mean1 == 0) {
    mean1 = 1
  }
  mean2 <- mean(as.numeric(num_list[cond_2]))
  if (mean2 == 0) {
    mean2 = 1
  }
  sd1 <- sd(as.numeric(num_list[cond_1]))
  sd2 <- sd(as.numeric(num_list[cond_2]))
  sd1 <- min(sd1, 0.2*abs(mean1))
  sd2 <- min(sd2, 0.2*abs(mean2))
  s2nvalue <- (mean1-mean2)/(sd1+sd2)
  return(s2nvalue)
}

protein_DGE$s2n <- apply(protein_quant,1,s2n,cond_1 = c(8:12),cond_2 = c(3:7))
# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_H_v5p2.rdata")
pathways <- Mm.H

# GSEA on Hallmark gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
  gseadata <- mirs.peaks[[i]]$target_Entrez_ID
  gseadata <- gseadata[!is.na(gseadata)]
  ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
  names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
  ranks <- sort(ranks, decreasing = T)
  fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
  fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
  suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
  fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
  if (dim(fgseaRes)[1] > 0) {
    GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_Hallmark_mergedpeaks.csv", sep = ""))
    for (n in 1:length(fgseaRes$pathway)) {
      p <- plotEnrichment(pathways[[fgseaRes$pathway[n]]], ranks) + labs(title = paste(mirna[i],"-",fgseaRes$pathway[n], sep = ""))
      grid.newpage()
      print(p)
    }
    grid.newpage()
    plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5) 
  }
}

pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_Hallmark_merged_peaks.pdf',
    width = 10,
    height = 6)
for (i in 1:20) {
  gseadata <- mirs.peaks[[i]]$target_Entrez_ID
  gseadata <- gseadata[!is.na(gseadata)]
  ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
  names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
  ranks <- sort(ranks, decreasing = T)
  fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
  fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
  suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
  fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
  if (dim(fgseaRes)[1] > 0) {
    GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_Hallmark_mergedpeaks.csv", sep = ""))
    for (n in 1:length(fgseaRes$pathway)) {
      p <- plotEnrichment(pathways[[fgseaRes$pathway[n]]], ranks) + labs(title = paste(mirna[i],"-",fgseaRes$pathway[n], sep = ""))
      grid.newpage()
      print(p)
    }
    grid.newpage()
    plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5) 
  }
}
dev.off()
## quartz_off_screen 
##                 2
names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.Hallmark.result.rds")

C2-curated gene sets

# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_c2_v5p2.rdata")
pathways <- Mm.c2

# GSEA on C2 gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
  gseadata <- mirs.peaks[[i]]$target_Entrez_ID
  gseadata <- gseadata[!is.na(gseadata)]
  ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
  names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
  ranks <- sort(ranks, decreasing = T)
  fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
  fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
  suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
  fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
  if (dim(fgseaRes)[1] > 0) {
    GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C2_mergedpeaks.csv", sep = ""))
    topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
    topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
    topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
    print(mirna[i])
    grid.newpage()
    plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5) 
  }
}
## [1] "let-7-5p/miR-98-5p"

## [1] "miR-15-5p/16-5p/195-5p/322-5p/497-5p"

## [1] "miR-29-3p"

## [1] "miR-200-3p/429-3p"

## [1] "miR-196-5p"

## [1] "miR-26-5p"

## [1] "miR-19-3p"

## [1] "miR-30-5p/384-5p"

## [1] "miR-484"

pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_C2_merged_peaks.pdf',
    width = 10,
    height = 6)
for (i in 1:20) {
  gseadata <- mirs.peaks[[i]]$target_Entrez_ID
  gseadata <- gseadata[!is.na(gseadata)]
  ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
  names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
  ranks <- sort(ranks, decreasing = T)
  fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
  fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
  suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
  fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
  if (dim(fgseaRes)[1] > 0) {
    GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C2_mergedpeaks.csv", sep = ""))
    topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
    topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
    topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
    print(mirna[i])
    grid.newpage()
    plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5) 
  }
}
## [1] "let-7-5p/miR-98-5p"
## [1] "miR-15-5p/16-5p/195-5p/322-5p/497-5p"
## [1] "miR-29-3p"
## [1] "miR-200-3p/429-3p"
## [1] "miR-196-5p"
## [1] "miR-26-5p"
## [1] "miR-19-3p"
## [1] "miR-30-5p/384-5p"
## [1] "miR-484"
dev.off()
## quartz_off_screen 
##                 2
names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.C2.result.rds")

C3-motif gene sets

# load the geneset data
load("~/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_09282019/Analysis/CLIP_Analysis/Data Visualization/GSEA Analysis/mouse_c3_v5p2.rdata")
pathways <- Mm.c3

# GSEA on C3 gene set for the top 20 miRNAs
GSEA.list <- list()
mirname <- c()
plist <- c()
for (i in 1:20) {
  gseadata <- mirs.peaks[[i]]$target_Entrez_ID
  gseadata <- gseadata[!is.na(gseadata)]
  ranks <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$s2n
  names(ranks) <- protein_DGE[protein_DGE$Entrez_ID %in% gseadata, ]$Entrez_ID
  ranks <- sort(ranks, decreasing = T)
  fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)
  fgseaRes <- fgseaRes[fgseaRes$padj < 0.1,]
  suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
  fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
  if (dim(fgseaRes)[1] > 0) {
    GSEA.list <- c(GSEA.list, list(as.data.frame(fgseaRes)))
    mirname <- c(mirname,mirna[i])
    write.csv(as.data.frame(fgseaRes), paste("GSEA Analysis/", str_remove_all(mirna[i], "/"), "_C3_mergedpeaks.csv", sep = ""))
    topPathwaysUp <- fgseaRes[ES > 0][head(order(padj), n=10), pathway]
    topPathwaysDown <- fgseaRes[ES < 0][head(order(padj), n=10), pathway]
    topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
    print(mirna[i])
    grid.newpage()
    plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5) 
  }
}
## [1] "let-7-5p/miR-98-5p"

## [1] "miR-29-3p"

## [1] "miR-200-3p/429-3p"

## [1] "miR-17-5p/20-5p/93-5p/106-5p"

## [1] "miR-24-3p"

## [1] "miR-27-3p"

## [1] "miR-103-3p/107-3p"

## [1] "miR-196-5p"

## [1] "miR-19-3p"

## [1] "miR-141-3p"

## [1] "miR-22-3p"

names(GSEA.list) <- mirname
saveRDS(GSEA.list, "Datafiles/merged.peak.GSEA.C3.result.rds")

Plot a heatmap for GSEA analysis done on all 2 genesets

# compile a matrix containing all genesets that had enrichment for all 20 miRNAs and their enrichment score
hall.GSEA <- readRDS("Datafiles/merged.peak.GSEA.Hallmark.result.rds")
c2.GSEA <- readRDS("Datafiles/merged.peak.GSEA.C2.result.rds")
#c3.GSEA <- readRDS("Datafiles/merged.peak.GSEA.C3.result.rds")

gsea.matrix <- data.frame(c(1:20))
gsea.matrix <- t(gsea.matrix)
colnames(gsea.matrix) <- mirna[1:20]
gsea.matrix[1,] <- NA
for (i in 1:length(hall.GSEA)) {
  for (n in 1:length(hall.GSEA[[i]]$pathway)) {
    if (length(grep(hall.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
    col.index <- grep(names(hall.GSEA[i]), colnames(gsea.matrix))
    row.index <- grep(hall.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
    gsea.matrix[row.index, col.index] <- hall.GSEA[[i]]$ES[n]
    }
  else {
    gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
    rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- hall.GSEA[[i]]$pathway[n]
    col.index <- grep(names(hall.GSEA[i]), colnames(gsea.matrix))
    gsea.matrix[dim(gsea.matrix)[1], col.index] <- hall.GSEA[[i]]$ES[n]
    }
  }
}

gsea.matrix <- gsea.matrix[-1,]

for (i in 1:length(c2.GSEA)) {
  for (n in 1:length(c2.GSEA[[i]]$pathway)) {
    if (length(grep(c2.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
    col.index <- grep(names(c2.GSEA[i]), colnames(gsea.matrix))
    row.index <- grep(c2.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
    gsea.matrix[row.index, col.index] <- c2.GSEA[[i]]$ES[n]
    }
  else {
    gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
    rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- c2.GSEA[[i]]$pathway[n]
    col.index <- grep(names(c2.GSEA[i]), colnames(gsea.matrix))
    gsea.matrix[dim(gsea.matrix)[1], col.index] <- c2.GSEA[[i]]$ES[n]
    }
  }
}

#for (i in 1:length(c3.GSEA)) {
#  for (n in 1:length(c3.GSEA[[i]]$pathway)) {
#    if (length(grep(c3.GSEA[[i]]$pathway[n], rownames(gsea.matrix))) > 0) {
#    col.index <- grep(names(c3.GSEA[i]), colnames(gsea.matrix))
#    row.index <- grep(c3.GSEA[[i]]$pathway[n], rownames(gsea.matrix))
#    gsea.matrix[row.index, col.index] <- c3.GSEA[[i]]$ES[n]
#    }
#  else {
#    gsea.matrix <- rbind(gsea.matrix, rep(NA,20))
#    rownames(gsea.matrix)[dim(gsea.matrix)[1]] <- c3.GSEA[[i]]$pathway[n]
#    col.index <- grep(names(c3.GSEA[i]), colnames(gsea.matrix))
#    gsea.matrix[dim(gsea.matrix)[1], col.index] <- c3.GSEA[[i]]$ES[n]
#    }
#  }
#}

gsea.matrix <- as.matrix(gsea.matrix)
my_palette <- colorRampPalette(c("blue", "white", "red"))(128)
na_color <- colorRampPalette("grey")

png('GSEA Analysis/GSEA_HEAP-CLIP_merged_peaks.png',
    width = 1200,
    height = 800,
    res = 100,
    pointsize = 8)
par(cex.main=1.8)
heatmap.2(gsea.matrix,
          main = "GSEA for 20 most active\nmiRNAs from merged HEAP-CLIP",
          density.info = "none",
          key = TRUE,
          lwid = c(1,7),
          lhei = c(1,7),
          col=my_palette,
          labRow = rownames(gsea.matrix),
          margins = c(20,25),
          symbreaks = TRUE,
          trace = "none",
          Rowv=FALSE, 
          Colv=FALSE,
          dendrogram = "none",
          na.color = "grey",
          cexCol = 1,
          na.rm = FALSE)
dev.off()
## quartz_off_screen 
##                 2
pdf('PDF_figure/KEGG_and_GSEA_analysis_merged/GSEA_HEAP-CLIP_merged_peaks.pdf',
    width = 24,
    height = 16)
par(cex.main=1.8)
heatmap.2(gsea.matrix,
          main = "GSEA for 20 most active\nmiRNAs from merged HEAP-CLIP",
          density.info = "none",
          key = TRUE,
          lwid = c(1,7),
          lhei = c(1,7),
          col=my_palette,
          labRow = rownames(gsea.matrix),
          margins = c(20,25),
          symbreaks = TRUE,
          trace = "none",
          Rowv=FALSE, 
          Colv=FALSE,
          dendrogram = "none",
          na.color = "grey",
          cexCol = 1,
          na.rm = FALSE)
dev.off()
## quartz_off_screen 
##                 2
GSEA

GSEA